Spin Structure Factor of the Frustrated Quantum Magnet CS2CUCI4 
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The ground state properties and neutron structure factor for the two-dimensional antiferromag- 
net on the triangular lattice, with uni-directional anisotropy in the nearest-neighbor exchange cou- 
plings and a weak Dzyaloshinskii-Moriya (DM) interaction, are studied. This Hamiltonian has been 
used to interpret neutron scattering measurements on the spin 1 /2 spiral spin-density- wave system, 
CS2CUCI4, [R. Coldea, et al, Phys. Rev. B 68, 134424 (2003)]. Calculations are performed using a 
1/S expansion, taking into account interactions between spin- waves. The ground state energy, the 
shift of the ordering wave- vector, Q, and the local magnetization are all calculated to order 1/S 2 . 
The neutron structure factor, obtained using anharmonic spin- wave Green's functions to order 1/S, 
is shown to be in reasonable agreement with published neutron data, provided that slightly different 
parameters are used for the exchange and DM interactions than those inferred from measurements 
in high magnetic field. 
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I. INTRODUCTION 

There is enormous interest in the condensed mat- 
ter community in the possibility of observing fractional- 
ized quasiparticles in two or three dimensional electronic 
systemsi^ Fractionalization is the rule, rather than the 
exception, in one dimensional metals where electronic 
spin and charge propagate independently. 4 Fractionally 
charged quasiparticles have also been found in the two di- 
mensional electron gas in high magnetic fields where the 
fractional quantum Hall effect is observed^ On the the- 
oretical side, fractionalization has recently been argued 
to exist in frustrated two-dimensional quantum model 
systems*^ although no corresponding experimental sys- 
tems have yet been confirmed. A set of long- wavelength 
theories have been developed that describe the properties 
of putative fractionalized magnets AsLiii Subsequently, 
the theoretical modeling of possible fractionalized phases 
was extended to a wide variety of strongly correlated elec- 



tron systems 
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The discovery of fractionalization in 



real two or three dimensional systems in zero magnetic 
field would provide a striking example of the emergence 
of new kinds of particles from the collective behavior of 
strongly interacting electrons^ 

As a result of this intense interest, the recent claim of 
the observation of spinons, neutral, spin 1/2 quasipar- 
ticles, in the insulating antiferromagnet CS2CUCI4 has 
attracted considerable attentionii 4 * 1 ^ CS2CUCI4 consists 
of nearly uncoupled layers of coupled chains. The chains 
are staggered, so that each spin interacts equally with 
two spins on each of two neighboring chains, with an an- 
tiferromagnetic coupling, J', about 1/3 the strength of 
the intrachain antiferromagnetic nearest neighbor cou- 
pling, J. The layer thus resembles a triangular lattice, 
with strong interactions along one direction and weaker 
interactions in the transverse directions. These intralayer 



couplings, J and J' assisted by a much weaker interlayer 
coupling J", stabilize an incommensurate spin density 
wave (SDW) ground state with a wave vector oriented 
along the chain direction. A weak Dzyaloshinskii-Moriya 
(DM) interaction D is believed to orient the spins in the 
2D layer. 

Coldea and coworkers 1 ^ 1 ^ used neutron scattering to 
study spin-wave excitations from this ordered ground 
state. While Coldea et al. did observe excitations sim- 
ilar to what is expected for spin-waves, they also saw 
broad features which they interpreted as excitations of 
pairs of spinons. Their idea is that the higher energy 
spin excitations resemble excitations from a nearby-lying 
spin-liquid state that supports fractionalized excitations. 
The low-dimensional nature of the system together with 
the low spin are argued to place this system close to 
a quantum critical point separating the ordered SDW 
state from a 2D spin liquid state. This conjecture has 
stimulated a number of theoretical studies of possible 
2D spin-liquid states that could account for the observed 
behavior, 16 i 17 . 18 i 19 i 20 i 21 . 22 

In this paper, we pursue a different interpretation of 
the neutron data, namely that low spin and quasi-one- 
dimensionality, along with non-collinear SDW order, all 
give rise to substantial anharmonic interactions that cou- 
ple one- and two-spinwave excitations from the ordered 
state. The coupling to two spinwave states also yields 
broad spectra as seen in the neutron data. A similar con- 
clusion has been reached by Veillette et al.^ 3 . who also 
noted that the low momentum resolution of the neutron 
data is another significant factor in the observed breadth 
of the spectra. We find that anharmonic effects lead to 
significant broadening of the neutron spectra but are sen- 
sitive to anisotropies, such as the DM interaction, which 
suppress quantum fluctuations. We show in particular 
that, if such anisotropies are weak, then quantum flue- 
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tuations, described in terms of a 1/S expansion about 
the mean field SDW state, reduce the local moment and 
lead to substantial renormalizations and broadening of 
the excitation spectra. 

The rest of this paper is organized as follows. In Sec. 
II we describe the Hamiltonian and the SDW state of 
CS2CUCI4. In Sec. Ill we show how anharmonic effects 
modify the ground state energy and the wave vector Q of 
the SDW to order 1/S 2 , and we also review the calcula- 
tion of the anharmonic one-spin-wave Green's functions 
to order 1/S. In Section IV we calculate the sublattice 
magnetization M to order 1/S 2 as a function of the ratio 
J'/ J and strength of DM interaction D/ J, showing that 
the renormalized value of M depends sensitively on the 
ratio D/J for small D/J. In Sec. V we examine expres- 
sions for the anharmonic spin- wave energies and damping 
to order 1/S. For D = we study how the Goldstone 
mode at Q is preserved and show how the preservation 
of this zero energy mode can be used to define a set of 
renormalized coupling constants. Section VI reviews the 
calculation of the neutron structure factor and presents a 
detailed comparison to the inelastic data of Coldea et al. 
for specific wave vectors where two-magnon scattering is 
important. In Sec. VII we discuss the implications of 
our theoretical calculations in interpreting existing data 
on CS2CUCI4. We also make suggestions for possible new 
experiments. 

II. SPIN DENSITY WAVE STATE FOR Cs 2 CuCl 4 

The spin density wave state of CS2CUCI4 can be de- 
scribed by the simple, model Hamiltonian for a set of 
decoupled layers: 

H = 2^ ly^R • Sr +< 5 1+< 5 2 + J'Sr • (Sr+^j + Sr + 5 2 ) 
R 

- (-1)"D.S r x(Sr +5i +Sr +52 )], (I) 

where J is the nearest neighbor coupling constant be- 
tween S = 1/2 Cu 2+ spins along chains in one direction, 
which we take to be the x-direction in a triangular lat- 
tice, while J', is the coupling constant along the other 
two principal directions in each layer as illustrated in 
Fig-H Both J and J' are antiferromagnetic, and J > J'. 
The last term in Eq. (QJ with D = (0, D, 0) describes the 
Dzyaloshinskii-Moriya interaction that alternates in sign 
between even and odd layers labeled by index n. 

Values of the coupling constants, determined from 
measurements in high magnetic fields^4 are: 

J = 0.374 meV, J' = 0.128 meV « J/3, (2) 

D = 0.02 meV, (3) 

leading to a small Neel temperature Tjv of less than IK 
in the presence of the interlayer coupling J" . The inter- 
layer coupling is sufficiently small (J" / J = 0.045) 24 that 




FIG. 1: Exchange couplings between the different sites of the 
triangular lattice within a single layer. 

we neglect it in our calculations performed at T = 0. 
The classical ground state for this Hamiltonian is a SDW 
whose wave vector Qo is along the strong-coupling direc- 
tion and has the value 

Qo = 27r(0.5 + e), (4) 

where experimentally, e = 0.030(2) is found for 
Cs2CuCi4r4 so that Qo is close to n. 

This means that individual chains are ordered approx- 
imately antiferromagnetically, while the spins in adjacent 
chains, which are slipped by half a lattice spacing along 
the chain direction, are approximately orthogonal. 

To further illustrate this point, consider the calculation 
of the classical ground state energy of the SDW state, 
that can be easily determined by introducing a local ref- 
erence frame, so that at every site the averaged spin is 
directed along the ( axis: 

Sr = 5"rCos(Q • R) - S R sin(Q • R), 
Sr. = 5 , R sin(Q • R) + S R cos(Q • R), 

S R = S^. (5) 

In Eqs. l[5)l. Q is the wave vector of the spiral structure, 
which at the classical level we write as (Qo,0,0) where 
Qo changes sign from layer to layer. Qo is determined 
from minimization of the classical ground state energy 
with S^ = 5 R = for all R. The classical ground state 
energy (per spin) is given by 

S 2 E£\Q) = S 2 J^ (6 ) 
Jq = Jq-iD q , (7) 

where 

k x \fik 

Jk = J cos k x + 2 J' cos cos — > (8) 

and 

At = -2tD sin — cos — (9) 
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are the Fourier transforms of the exchange and DM in- 
teractions of the Hamiltonian respectively. 
For D = 0, e in Eq. is given by 



e = - arcsin ( ^— I = 0.0547, 

7T V 2J , 



(10) 



while for D = 0.02 meV, the value of Qq that mini- 
mizes the classical ground state energy corresponds to 
e = 0.0533 for the experimental values of J and J'. The 
two terms in Eq. J5J), represent the intrachain and inter- 
chain energies respectively. Substituting the actual value 
of Qo for a given J'/ J and D = 0, we find that 



S 2 E%\q ) = -JS 2 



1 



2J 2 



(11) 



Thus the contribution to the ground state energy from 
the interchain coupling is of order (J') 2 / J 2 compared to 
the intrachain contribution. Quantum fluctuations re- 
duce the effectiveness of interchain coupling further by 
renormalizing the SDW wave vector to a value even closer 
to 7r (corresponding to a smaller effective J'/ J). At the 
same time, quantum fluctuations lower the total energy. 
Thus if one were to model the effects of quantum fluc- 
tuations in terms of effective coupling constants, J and 
J', the result would be J' /J < J' / J and J > J. As is 
shown in the next section and in Appendix A, the value 
of e obtained at the classical level, is in fact considerably 
renormalized by quantum fluctuations. 



III. 1/S EXPANSION 

To address the physics determined by quantum fluc- 
tuations, we employ the well-known Holstein-Primakoff 
transformation for the spin operators^ 



■Sr. — "Sr. 



Sr — Sr 



R 



2Sa 



R 



1 



1 - ^ a R«R 



45' 



OrOr 



(12) 



(13) 



S R = S ~ a R a R- 



(14) 



These are written in a form with the square root ex- 
panded to first order in 1/(25*). The magnon operators 
a R,i a R describe excitations around the spiral ground 
state and obey Bose statistics. Substituting the trans- 
formations 112ll - 114|) into Eq. we obtain the Hamil- 
tonian for interacting magnons, 

n = s 2 e < §\q) + (h (2) + n {3) + h {4) ), (15) 

where in the Fourier transformed representation 



H (2) = 25^T ( A k al a k - -^(4a f _ k + a k a_ k ) ) . (16) 



2 

The functions A k and are expressed through J k as 
1 



Q+k 



J, 



Q-kJ ■ 



(17) 



(18) 



The interactions between magnons are described by 
the last two terms Ti^ and Ti^ . The three-magnon 
term can be written in the form 

n (3) = iJJ-^id + C 2 )(4a 2 ai 

1,2,3 

- a{a f 2 a 3 )S 1+2 .3 (19) 



that contains the vertex 

Ck = (j£ +k - j£_ k ) , (20) 

which is antisymmetric with respect to the transforma- 
tion k — > — k. For brevity, in longer expressions we use 
the convention that 1 = ki, etc. 

The four-magnon term can be conveniently written as: 



nii) = J2 { t( Al - 3 + Al - A + A2 - 3 + A2 -^ - ( Bl - 3 + Sl ~ 4 + B2 - 3 + B2 -^ - ( Al + A 2+A 3 + At) 

1,2,3,4 

a\a 2 a 3 a 4 5 1+2 , 3+4 + ~ (Bi + B 2 + B 3 ) {a\a\a\a A + a\a 3 a 2 a^j <5i +2 +3,4 j , 



(21) 



and contains only symmetric vertices. It is important to 
emphasize that TL^ comes from the coupling between 
the operators 5 R and 5 R , which arises for non-collinear 
ordered states. This term plays a crucial role in calcula- 



tions of the structure factor as well as in renormalization 
of the energy spectrum. 

To proceed, one needs to diagonalize the quadratic part 
of the Hamiltonian Eq. (|16fl by means of the Bogoliubov 
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transformation: 

a k = M kC k + U k ct k , A k = w k cl k + v k c k , (22) 



u k = 



:k D / A k - e k 



(23) 



2e k V 2e k 

The energy spectrum in the above equations is given by 



?k 



(Jk-Jg) 



(-/, 



T i tT 



Q+k 



Q-k 



)/2-^ 



(24) 



For D = 0, the magnon spectrum has zeros at k = 
and k = ±Q in the two-dimensional Brillouin zone, while 
a non-zero D leads to the appearance of a finite gap of 
order V £> J around the Q-points. After diagonalization, 
7i( 2 ) takes the form: 



iV 



1 



k 



The contribution Se£\q) gives the leading 1/(25) 
correction to the classical ground state energy Eq. JJjJ. 
Minimizing the sum S 2 E { °\Q) + SEg\Q), we find the 
quantum correction to the classical value Qq. To order 
1/5 the shifted wave vector is Q = (Q, 0, 0) 



(25) 



(26) 



Q = Qo + 



25 



AQ 



(i) - 



5Q 2 



1^(4 + s k ) 3J^ +k 
2V ^ 

k 



£k 



9Q 



In the absence of the DM interaction, 
AQ (1) /2tt = -0.0324, 

while for D = 0.02 meV, 

AQ (1) /2tt = -0.0228. 



(27) 
(.28) 

'o 

(29) 
(30) 



We see that inclusion of the Dzyaloshinskii-Moriya inter- 
action suppresses the renormalization of Q towards it. 

It is possible to go further and also calculate the cor- 
rection to the ground state energy and the ordering wave 
vector Q that is of order l/(25) 2 . The details of those 
calculations are presented in Appendix A. The results in- 
dicate that quantum corrections computed order by or- 
der in 1/(25) are significant for the values of J and J' 
given by Eqs. (J2J. To further explore this issue, one 
can compare the classical ground state energy Eq. JBJ to 
its leading renormalization due to quantum fluctuations 
SE$p(Qo). We find, for D = 



S 2 E { °\q )/J = -0.265, SE£\Qo)/J = -0.157,(31) 



and, for D = 0.02 meV, 



S' z E { °\q )/J = -0.291, SE£'(Q Q )/J = -0.138.(32) 

We see that the leading quantum correction lowers the 
ground state energy by about 50%, giving reasonable 
agreement with a recent experimental determination of 
this quantity, based on susceptibility measurements, 26 
which yields a ground state energy slightly above —0.5 J. 
We show in Appendix A that for D = the next order 
in 1/(25) correction reduces the ground state energy by 
a further 10 percent below the sum of the values given in 
Eq. (|31|l . bringing the total within a few percent of the 
experimental value. 

To calculate physical observables, we need the Green's 
functions for the magnon operators ak. They can be 
conveniently written in the form of a 2 x 2 matrix, indi- 
cating the presence of normal and anomalous parts. At 
zero temperature, the definition reads 



(!)/ 



iG(k,«) 



dte 1 



T 



-kW 



4(0)a_ k (0) 



(33) 



where T is the time ordering operator, and the average 
is taken with respect to the ground state. The inverse of 
the matrix for the bare Green's functions has the form 29 



«s™-'<M>-("-££ + *^jjSL + A (34) 



corresponding to the Green's function 



G(°)(k,uA = - ^ (" + 2 ^-iry 255 y 

v ' (to - + tr))(uj + low - ir/) \ 25B k -u + 2SA k - vq ' 
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where lo-^ = 2Se k . The self-energy S(k, oS) determining the exact Green's function is also a 2 x 2 matrix, that can be 
parametrized as 



f( , A _ /x(k,«) + y(k, w ) 



and which satisfies the Dyson equation 

G _1 (k,w) = G^-Hk.w) - E(k,w). (37) 

The self-energy to order 1 / (25) consists of two parts 

£(k,w) = £( 4 >(k) + £( 3 >(k,w). (38) 

The contribution E^ 4 ) (k) is frequency independent and 
can be calculated simply by the Hartree-Fock decoupling 
of the quartic term in the Hamiltonian Eq. Ij21|l . 



JfW(k) 



25 



E- 



-Bic + B„ B, 



(A k _ q - B k - q - Aq - A k ) A q 



(39) 



Z(k,w) 
X(k, w) - F(k, w) 



(36) 



r 



y( 4 )(k) 



25 

AT 



E- 



k-q — q ^ • 'q 



.4 



— B k -q — A q — -A k -B, 



Z< 4 >(k) = 0. 



(40) 



(41) 



E( 3 )(k ,u>) is the contribution arising from the three- 
magnon interactions, Ti^ . Its non-zero imaginary part 
determines the magnon damping to first order in 1/(25). 
The corresponding components are most easily calculated 
by transforming Eq. I|19l) to the Ck-operator basis using 
the Bogoliubov transformations, Eqs. I|22|). 



X (3) (k,w) 



5 
167V 



^{[a> 1 (q,k-q)] 2 + [<f>2(q,k-q)] 2 } 



+ 



Wq + W k -q — W - IT\ Wq + W k -q + LO - ITj 



Z^(k,u;) = --^^{^(q.k-qjf-^q.k-q)] 2 } 



1 



1 



Wq + ^k-q - W — IT] LOq + Wk-q + UJ — 11] 



S 



16iV 

$i(q,k-q) 
$ 2 (q,k-q) 



1 



Wq + Wk-q - LU — IT] LUq + Wk-q + W - 27? 



(C q + C k -q) (tlq + Uq) (Ufc-q + «k-q) - 2C k (Uq«k-q + ^q^k-q) , 
Cq (lt q + Wq) (ttk-q ~ «k-q) + Cfc-q («k-q + «k-q) ("q ~ «q) ■ 



(42) 



(43) 



IV. SUBLATTICE MAGNETIZATION 



In this section we consider the staggered magnetiza- 
tion for a range of values of J' / J and D/J and show that 
quantum fluctuations lead to a considerable renormaliza- 
tion of the average value of the local moment. From the 
operator definition Eq. (|14f) we can write the staggered 
magnetization as 



M = S + 



N 4-^ 



2^1 



G u (k,£ i ;)e" 



(44) 



where the exact Green's function Gn(k, oj) is the corre- 
sponding element of the full matrix G(k, u>) and 0_ is 
a negative infinitesimal. To lowest order in 1/(25), one 

can use G^(k, w) of Eq. I|35|) . to find the first order cor- 
rection to the 5=1/2 value 



(45) 



For the next order correction in 1/(25), there are two 

(2) 

contributions. The first one, denoted as A.M} , arises 
from renormalization by quantum fluctuations of the or- 
dering wave vector Qq. This correction is obtainable by 
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substituting the renormalized value Q given by Eq. I(27() The second contribution is determined by the self-energy 
into the formula above and expanding using Eq. 145(1 up corrections of order 1/(25) to the Green's function itself 
to first order in 1/(25*). and is given by 



AM 



(2) 



25 N 



lei 



dQ 



(46) 



Qo 



AM, 



-(2) _ 1 f 



dio 



-iu>0_ 



AGn(k,o;), 



(47) 



where 



AGn(k, ui) = (G(°)(k , w)£(k, oj)G^ (k, u) 
[X(k, u>) + Y (k, w )] [gJ? (k, w)] 2 + 2Z(k, W )G£? (k, wJG^ (k, «) + [X(k, u) - Y(k, w )] [g$ (k, w) 
This correction, as a result of integration over u>, will contain two parts, so that 

AM : (2) = AM' 2 ' + AMf. 



(48) 



(49) 



One part comes from taking the residues at the poles of self-energies, while the other results from the double-poles of 
the products of two Green's functions. Using Eqs. ((391) - ((42(1 . we obtain the first contribution 



AM? } 



S \ ^ 
~8N 2-*> 



[*i(q,k-q) 



G^(k, fi k , q ) + Gi l 2 ; (k,O k , q 



(0), 



k,q 



G<° } (k, n k , q ) - G[° 2 ] (k, O k , q ) - (q, k - q)* 2 (q, k - q 



2 , [* 2 (q,k-q)r 



G£ 3 (k,ft k , q ) - Gg)(k,n k , q 



(50) 



k,q 



-UJq — W k _ q . 



(51) 



The contribution arising from the residues at the double-poles of the products of two Green's functions in Eqs. 
(147(1 - 1(4811 . is straightforward to calculate as well. After some algebra, we find that 

AM< 2) = -i- ^ || | [X(k, + Z(k, c*)] [A k + B k ] - [X(k, - Z(k, c*)] [A k - S k ] | 
, 1 v^(4- £k) 



JV 



44 



[X'(k, wk) + Z'(k, wt)] [A k + B k ] - [X'(k, cu k ) - Z'(k, wk)] [A k - S k ] - 2Y'(k, Wk )e k k (52) 



where w k = 25e k , and 



\ duj 

r 



(53) 



(2) 

It is clear from the expressions above that both AMj 

(2) 

and AMjj contain the overall prefactor 1/(25). 

In Fig. |2 we plot the total sublattice magnetization 



ratios of each correction to the classical magnetization 
5 = 1/2, —AM/5, are also plotted, as explained in the 
caption of Fig. One can see that the renormalization 
of the local magnetization by quantum fluctuations to 
second order in 1/(25) is considerable for all values of D 
between and 0.025 meV, although increasing D clearly 
as a function of Dzyaloshinskii-Moriya interaction D, for suppresses fluctuations making the system more classical. 



M = S + AM (1) + AMj (2) + AMj (2) (54) 

:tion of Dzyaloshinskii-Moriya intera 
coupling constants J and J' as given by Eqs. 0. The The n S ure also shows a calculated sublattice magnetiza- 
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-0.05 



0.000 



0.005 



0.010 0.015 

D (meV) 



0.020 



FIG. 2: The total value of sublattice magnetization M, Eq. 
I54L calculated up to the second order in 1/(25) as a function 
of Dzyaloshinskii-Moriya interaction D, is shown by the solid 
line. The dashed, dotted and dashed-dotted lines are the 
results for the quantum corrections -AM' 1 ', -AM] 2 ' and 
— AMjj 2 ', respectively, divided by S = 1/2. 



tion that is negative for D between and D 



0.003 



meV, signaling the melting of the assumed spiral ordered 
ground state by quantum fluctuations when D is too 
small. 

Experimentally, a value of —AM/5 s» 0.25 was in- 
ferred by Coldea and co-workers^I from elastic neutron 
scattering, and this value was used in the analysis of the 
most recent neutron data^ Such a small renormaliza- 
tion of the local moment is incompatible with the the- 
ory presented here. It is also unlikely to be explained 
by the addition of the very small expected interlayer ex- 
change coupling J" I J = 0.045ii2i2i It is indeed puzzling 
that such a small moment reduction, which is less than 
for the spin 1/2 Heisenberg antiferromagnet on a two- 
dimensional square lattice, would occur in this system 
which is both frustrated and close to the one-dimensional 
limit. 



V. ENERGY SPECTRUM 

The rcnormalized magnon energy spectrum ci) k — 25e k 
can be found from the poles of the exact Green's function, 
determined by2& 



He{detG _1 (k,wi t )} = Re|det G< > _1 (k,Wk) - S(k, <D k ) } = 0, 
which leads to the self-consistency equation: 
Wk = 



(55) 



{{2S{A k + B k ) + Re [X(k, £ k ) - Z(k, £ k )]} {25(,4 k - B k ) + Re [X(k, cD k ) + Z(k, Co k )]} 
- {Im [X(k, u> k ) + Z(k, £>k)]} {Im [X(k, uj k ) - Z(k, £ k )]} + [ImY(k, erf} 17 ' + R C y (k, cD k ). (56) 

Near these poles, if the imaginary part of the self-energy is small, the shape of the one-magnon spectrum is approxi- 
mately Lorentzian and has the formSi 



detG 1 (k, oj) ~ — 2w k z k (w — w k ) 



ilm 



detG _1 (k,w k ) 



(57) 



r 



where 



OU) 



detG _1 (k,a;) 



(58) 



The imaginary part of Eq. (|57J) determines the inverse 
life-time of the quasiparticles. To calculate the half-width 
at half-maximum of the peak, we use the parameter 



ImtdetG-^k.Wk)] 



2cj k z k 



(59) 



Equation (|56|l . although formally exact, implicitly in- 
cludes contributions of all orders in 1/5. Furthermore, 



quantities such as the terms 25(A k ± _B k ) that are of 
the leading order in 5, should be corrected for the first 
order shift of the ordering wave vector from Qo to Q, 
given by Eqs. H27 jH28fl . Although the self-energies are 
only calculated to order 1/5, taking the square root of a 
sum of squares of self-energies evaluated at renormalized 
frequencies, effectively mixes in higher order corrections. 
Also, particularly for the case D = Q, there is cancellation 
among terms of order 1/5, leaving a result, at the shifted 
Q- vector, which is the square root of residual 1/5 2 terms, 
which may be negative. 

A simple way to preserve physical behavior at the 
Goldstone wave vector Q for D = 0, is to evaluate the 
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self-energies using renormalized coupling constants, for 
which the classical ordering wave vector is Q, that is 



D = 0.02 meV 



D = meV 



, J' 

7r + 2 arcsm | — = 
. 2j 



(60) 



The effect of this renormalization, which determines the 
ratio J' I J and represents a higher order in 1/5 correction 
to the self-energies, is that the self-energies will vanish at 
the renormalized Q. A further renormalization of J by 
itself is defined by the condition that the coefficient of 
the linear term in |k — Q| inside the square root vanishes 
for k = Q, so that the energy dispersion relation around 
this point is linear in |k — Q|. The procedure of expand- 
ing the self-energies around Q is tedious but completely 
equivalent to that described elsewhere pSLSi Omitting the 
intermediate steps, we obtain the condition 



<9k 



1 



k=Q 



257V 



E 



(A p + B p ) _ dJp+Q 
Ep dp 



(61) 



where the bars indicate that the renormalized values J 
and J' should be used in the right hand side of this equa- 
tion. One can easily check that, within this approach, we 
retain consistency up to lowest order in 1/25. Indeed, 
expanding the derivative dJ^/dh around Qo, we recover 
Eq. (|28|l for the correction to the ordering wave vector. 
Using Eqs. ©, we find from Eqs. fffl ty -jST ) that 



J w 0.472 meV, J' w 0.066 meV. 



(62) 



We see that in this approach, J' is renormalized down 
from its bare value, while J increases, as discussed in 
Sec. II. 

In spite of the fact that we expect the system to be 
unstable for D = and J'/ J — 1/3, it is of pedagogical 
interest to examine the effects of anharmonicity and the 
shift of the Goldstone mode for this case. The right hand 
panel of Fig. shows the results for the renormalized 
energy spectrum and its comparison to the linear spin- 
wave dispersion along the (1, 0) direction for D = 0. The 
inverse lifetime of quasiparticles Eq. (|59|l is also shown. 
The excitations are seen to be heavily overdamped close 
to the (10) point. Also kinks are seen at various places in 
the renormalized curves, apparently at points where one- 
to-two magnon decay channels turn on. Although not 
shown in Fig. |3| we also find that the band width of the 
dispersion along (1,0) — * (1, 1) is renormalized downward 
by about a factor of 3 due to anharmonic effects. 

In the left hand panel of Fig. we plot the renor- 
malized energy spectrum for D = 0.02 meV^ In the 
presence of this DM interaction, it suffices to use the 
bare values for the coupling constants in expressions for 
self-energies, since the finite gap in the spectrum does 
not lead to any problematic behavior near Q when D is 
not too small. The effects of quantum fluctuations are 
very sensitive to the presence of a non-zero D. Increasing 
the value of D reduces the interaction-induced damping 
of quasiparticles, making them more stable. Quantum 



0.4 



0.3 

> 

sg. 0.2 

3 

o.i 



0.0 



bare . ~* 

a ^y \ 






r \ \ ren /» 

/ A \ / \ 
/ \ \ 


o 




1 \ 1 ' n 


< 

a ■ 

a. 

sl" 

IB.- 

'■X - 

' '.° 


\ \ \ 

v \ 1 X rcn \ 

A// k 1 


/ \ 1 1 
\ \ \ 1 
1 v / ' 

/ r » 1 <k 

l k 1 ' \ 


: /' \ 

*i — | — i — 'i i 1 ■■ ■ i — i — 1 


\ \i ' ' 

1 \ ■ \ \ I 1 bare' 

V V '°\ .;• 
/ •' •\V 
I ' ' \\ 1 
■If' 


I- 
s . 



0.4 



> 



0.2 . 



1.0 0.8 0.6 0.4 0.2 0.0 0.2 0.4 0.6 0.8 1.0 



FIG. 3: (Color Online) The two panels present the results 
for the bare and renormalized energy spectra for the cases 
of D = and D = 0.02 meV. Spin-wave energies are shown 
together with the inverse lifetime Fk- The path in the two- 
dimensional Brillouin zone sweeps from (00) point to (10). k x , 
k y are measured in units of 2n and 2ty/V3 respectively. 



fluctuations, in addition to shifting Q towards w, reduce 
the size of the gap at this point. However, the effect of 
suppression of the gap towards zero, is relatively small 
for D = 0.02 meV. 



VI. SCATTERING INTENSITY 

In this Section, the dynamical properties of an anti- 
ferromagnet described by the model of Eq. (JTJ are cal- 
culated within the framework of the 1/(25) expansion. 
We are specifically interested in answering the question 
of whether anharmonic spin waves can account for the 
main features seen in inelastic neutron scattering mea- 
surements on CS2CUCI4.— We begin with expressions for 
the structure factor components that enter the formula 
for the inelastic neutron scattering cross section. These 
results were obtained earlier by Ohyama and Shiba, 29 
and we simply quote them here. 

Neutron scattering spectra of spin excitations in mag- 
netic solids may be expressed in terms of the Fourier 
transformed real-time dynamical correlation function, 
i,k = (x,y,z) 



1 r c 
5 lfe (k,^) = - / 



k-R) 



R 



(63) 



For the spiral spin density wave state, S xx (k,uj) = 
S zz (k,ui). The complete expression for the inelastic, dif- 
ferential scattering cross-section, including polarization 
factors, is given bj^ 

7(k,w) = Px S xx (k,Lu)+ Py Syy(k,u;), (64) 
p x = 1 + cos 2 a k , p y = sin 2 a k , (65) 
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where is the angle between the scattering wave vector 
k and the axis perpendicular to the plane of the spins. 

We begin by calculating the corresponding time- 
ordered spin-spin correlation functions in the rotated co- 
ordinate system 



iF af3 (k,Lu) 



dte- luJt (TS\(0)S£(t)), (66) 



where a, f3 = (£, T], £) are the rotated coordinate axes 
given by Eqs. (JSJ. The dynamical structure factor is re- 
lated to the imaginary parts of these correlation functions 
in the following way: 



with, 



5 OT/ (k,w) 
S xx (k,uj) 
S xz {k,Lu) 



-imF' ? ' ? (k,w), 

7T 



= S zz {k,u) = — Im[e+(k + Q ! w) + e-(k-Q,w)] , 



-S zx {k,Lu) 



e ± (k,c)- i {F«(k, £ 



-Im [0+(k + Q, u) - 0- (k - Q, w)] , 



F«(k 



(67) 
(68) 
(69) 

(70) 



We treat the interactions between the magnons as a perturbation in the formally small parameter 1/(25), and use 
the standard 5-matrix expansion while calculating the averages of the magnon operators a k (i)^& As a result, 



i™(k,w) = !c£[Gii(k,w) + G 2 a(k,w)-2Gi 2 (k,w)], 
F«(k,w) = |4[Gii(k,w) + G2a(k,w) + 2G 12 (k,w)], 



(71) 



where the Green's function are determined from Eq. i|37|) with the self-energies given by Eqs. fi fy -QB ^ . and, to the 
relevant order in 1/(25), 



1 



45N 



k + UkVk) 



(72) 



The term that mixes the transverse and longitudinal fluctuations has the form: 

i [F*< (k, u) - F« (k, w)] = {PW (k, «) [G n (k, w) + G 22 (k, w) - 2G 12 (k, w)] 



-P< 2 )(k,^) [G u (k,w)-G22(k,w)]} 



with the functions pt 1 ' 2 ) (k, u>) defined as 



P«(k,w)= — ^<& (1) (q,k-q)( Uq t- k _ q + ^ k _ q ) 



4iV 
5 



1 



1 



p( 2 )(k,c) = — ^$( 2 )( q ,k- q )( UqWk _ q + ^ k _ q ) 



Wq + UJk-q ~ LU — IT] LUq + UJk-q + LU — IT] 
1 1 



Wq + l^k-q — LU — 17] LUq + W k -q + LU — IT] 



(73) 



(74) 



As before, a-> k = 25e k . The longitudinal correlations expanded in inverse powers of 25 can be expressed as a sum of 
two contributions 



P«(k, W )=P Q cc (k ) a;)+P 1 cc (k, W ) 



where 



P Q CC (k, LU) = - ^ 51 ("q^k-q + WqWk-q) 



1 



1 



+ ^ k -q — LU — it] LUq + CJ k _ q + LU — 17] 



(75) 
(76) 
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Jf(k,w) = — |(pW(k,w)) [G 11 (k^)+G 22 (k,c) + 2G 12 (k^)] + (p( 2 )(k^)) • 

[G u (k,w) + G 22 (k,a;) - 2Gi 2 (k,w)] + 2P«(k, w )p( 2 )(k, w) [G u (k,w) - G 22 (k,w)]} . (77) 



Note that the term Ff (k, w) is formally of order 1/(2S) 2 
compared to the leading one magnon contributions Eq. 
(I71|l . However, it is kept because it contains the self- 
energy corrections to Eq. (|7fi|> describing the two magnon 
continuum. 

Using the set of formulas above, we have calculated 
the scattering intensities as a function of frequency for 
various trajectories in k and lu space, corresponding to 
the measured spectral In general our results agree with 
those of Veillette et ali 2 ^ except for small differences due 
to their use of self-consistent energy denominators in the 
self-energies and our use of shifted Q values in the anhar- 
monic Green's functions. Veillette et al. restricted their 
attention to values of J, J' and D given by Eqs. 
In the discussion that follows, we will also consider the 
consequences of varying the D parameter, which both 
shifts certain modes and also controls the size of quan- 
tum fluctuations, as well as varying the overall energy 
scale, J. 

We focus our attention on two trajectories in the (k, uj) 
space whose parameterizations correspond to the scans G 
and J of Ref . Q 

G: k x = 0.5, k y = 1.53- 0.32w - O.lw 2 

J : k x = 0.47, k y = 1.0 - 0.45w. (78) 

In Eq. Ij78(l , k x and k y are written in the units of 2ir and 
27r/v3 respectively, while the energy is measured in meV. 
Since the momentum transfer lies in the x-z plane, p x = 
p y = 1 for both of those scans. It is important also to take 
into account the finite energy and momentum resolution 
of the scattered neutrons. The finite energy resolution 
Aw w 0.016 meV leads to a small Lorentzian broadening 
of the one-magnon peaks. However, the effect of smearing 
due to the finite momentum resolution Ak/2ir — 0.085 
is much more pronounced and considerably broadens the 
structure of the G and J scans. What is special about 
the G and J scans is that they both exhibit sharp low- 
energy modes with energies around 0.1 meV. As we shall 
see, the energies of these peaks are quite sensitive to the 
exact value of the DM interaction D. 

In Fig. 01 we plot the results of scattering intensity 
calculations for scan G, a data set which was analyzed 
in detail by Coldea et al^» In particular, it was claimed 
that the high-energy tail in this scan has too much inten- 
sity, by nearly an order of magnitude, to be explained by 
linear spin-wave theory. Our calculations were done us- 
ing the bare couplings J and J', and the Dzyaloshinskii- 
Moriya interaction was set equal to D = 0.02 meV, the 
value obtained from experiments in high magnetic field. 24 
The theory predicts three peaks, one corresponding to 
the principal mode and the other two, reflecting the 
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FIG. 4: (Color Online) Intensity of the scattered neutrons for 
scan G in Eq. (ZB) i i n the presence of Dzyaloshinskii-Moria 
interaction D — 0.02 meV. The energy and momentum res- 
olutions are taken to be Aui = 0.016 meV and Afc = 0.085 
respectively. The thick solid line is the total scattering inten- 
sity, Eq. 1648 . Other lines represent the various contributions 
appearing in Eqs. 167(1 - 1711 . 



presence of the shifted secondary spin waves at £Q±k- 
For the non-renormalized couplings J and J', Eq. I|24|) 
for the bare spectrum predicts that the principal peak 
should be at 0.22 meV while the secondary peaks are 
very close to each other and located approximately at en- 
ergy 0.28 meV. Inclusion of self-energies to order 1/(25*) 
renormalizes the position of the main peak down to 0.18 
meV. The principal one-magnon peak should lie below a 
two-magnon scattering continuum starting exactly at the 
location of the principal peak. However, the presence of 
self-energies in the Green's functions, leads to a shift in 
positions of their poles, and to a finite gap between the 
principal peak and the start of the two magnon scatter- 
ing continuum. This result is an artifact of truncating the 
perturbation expansion at a finite order, as discussed by 
Veillette et ali^ We also observe that the two secondary 
peaks are so close to each other that it is difficult to dis- 
tinguish them, once the finite frequency and momentum 
resolution is fully incorporated into the calculations. Ex- 
perimentally, the main peak is observed at a somewhat 
smaller energy of 0.107 meV, while the secondary ones 
are located around 0.25 meV. 14 Fig.0]also shows the sign 
and magnitude of the various contributions to the total 
intensity. It is seen that the two-magnon part, integrated 
over the full range of frequencies, carries a considerable 
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FIG. 5: (Color Online) Intensity of the scattered neutrons 
corresponding to scan G, for the value of Dzyaloshinskii-Moria 
interaction D — 0.01 meV. Upper panel: The total structure 
factor with the energy (momentum) resolution AE = 0.016 
meV (Ak = 0.085). Inset: Experimental data from Ref. Il4l 
Lower panel: Different components which contribute to the 
structure factor (see Eqs. 167H - 17H in the text). 

overall weight compared to the one-magnon contribution. 
At the same time, the term mixing the transverse and 
longitudinal fluctuations can have both signs, and being 
small, does not influence the intensity profile very much. 

Why is the two-magnon intensity in our calculation 
so much larger, relative to the principal peak than in 
the analysis of Coldea et al.? The reason is the factor Cy 
which renormalizes the out-of-plane scattering and which 
is basically due to the reduction of the local moment by 
quantum fluctuations. As noted above, theory predicts 
a much smaller local moment than the surprisingly large 
value obtained from experimental According to our cal- 
culations, the moment reduction increases with decreas- 
ing D, and the ordered moment vanishes for D < 0.003 
meV. 

In Fig. [51 we plot the theoretical result for scan G, but 
this time with a reduced DM interaction D — 0.01 meV. 
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FIG. 6: (Color Online) Intensity of the scattered neu- 
trons corresponding to scan J, Eq. 1781 . in the presence of 
Dzyaloshinskii-Moria interaction D — 0.02 meV. The differ- 
ent components are displayed as in Figs. 21 03 



The comparison of our calculations with the experimen- 
tal measurements reported in Ref. is displayed in the 
upper panel. We see that for this value of D, the position 
of the principal magnon peak is shifted towards smaller 
energies. The gap between the main peak and the energy 
where the continuum appears, is also decreased; there is 
further broadening of the principal peak and more inten- 
sity in the two-magnon continuum. It appears that the 
value D = 0.01 meV provides a better explanation of 
the form of the total scattering intensity for scan G than 
D = 0.02 meV. Comparison of the plots for a range of 
values of D, shows that the scattering intensity is very 
sensitive to the strength of the DM interaction. [Note 
that the DM interaction was not included in the analysis 
of Ref. El] Finally, we observe that even this smaller 
value of D = 0.01 meV does not account for all of the 
high energy scattering which seems to extend up to 1 
meV, suggesting additional correlations at high energy 
which are not contained in the theory. 

Next, we examine scan J which probes the vicinity of 
the ordering wave vector^ Inelastic scattering at this 
point should show a primary peak at the wave vector 
kj w Q and two secondary modes, corresponding to mo- 
mentum transfer kj ± Q. The primary mode, which is 
prominent in the theoretical spectrum, directly measures 
the strength of t he D M interaction, since its energy is 
proportional to V DJ. The experimental data for point 
J, shown only in the inset of Fig. 5G in the paper of 
Coldea and co-workers^» do not extend to low enough 
energy to display this peak. Nevertheless the position of 
the peak, at an energy of 0.10(1) meV, is listed in Ta- 
ble 1 of that paper. This value of the mode at the J 
point is noticeably smaller than anharmonic spin-wave 
theory would predict for D = 0.02 meV, suggesting that 
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the effective value of the DM interaction constant D is 
smaller than 0.02 meV. Since the theory, for D = 0.02 
meV, gives a peak at 0.15 meV, a fit to the measured 
peak at 0.1 meV would suggest a value of D in the range 
of 0.008 to 0.01 meV. Of the two secondary peaks at the 
J point, one has an energy close to zero, corresponding 
to the shifted zone-center mode, while the other lies at 
higher energy, just above 0.3 meV. 



VII. DISCUSSION AND CONCLUSIONS 

The main motivation for the work presented in this 
paper was to see if a relatively conventional anharmonic 
spin-wave calculation could account for neutron scatter- 
ing measurements on CS2CUCI4. Our conclusion is a 
qualified "yes". 

Without any calculations, one knows that the observed 
excitations are from an ordered state exhibiting magnetic 
Bragg peaks, and hence that the lowest-lying excitations 
must be spin waves. However, both linear and anhar- 
monic spin-wave theory suggest a substantial moment 
reduction for a broad range around the expected values 
of Dj J and J'/ J. Further experiments would be useful 
to verify the amount by which the moment is reduced, 
since the present result^ is inconsistent both with an- 
harmonic spin- waves and with close proximity to a spin- 
liquid state. 

At low energies we find that we can fit the sharp spin- 
wave features of the data, albeit with a somewhat smaller 
value of the DM interaction constant, D, than was ob- 
tained from high field measurements of a fully polarized 
state. We also find from theoretical calculations that, 
without a sufficiently large DM interaction, long range 
order is destroyed by quantum fluctuations. However, 
the value of D that best fits theory to experiment is well 
within the range in which order is stable, at least for 
T = 0. Additional experiments, including measurements 
of the renormalized D in zero field, might clarify this 
point. (We note in passing that a value of J about 30- 
40% larger than given by Eq. J2Jl, for fixed J' /J, would 
provide a much improved fit to the higher energy one- 
spin- wave features.) 

Our theory yields a considerable amount of continuous 
two-spin-wave scattering. However the amount of this 
scattering, although readily measurable, is somewhat less 
than is observed experimentally. Furthermore, calcula- 
tion of the relative size of one- and two-spin-wave fea- 
tures is complicated by their sensitivity to the size of D. 
That is, a smaller D suppresses one-spin-wave scattering 
and makes two-spin-wave continuous scattering relatively 
more conspicuous. Quantum fluctuations also broaden 
the one-spin- wave features via lifetime effects. However, 
it is difficult to assess, from the data alone, how much of 
the observed broadening of the principal peaks is due to 
lifetimes and how much is due to low momentum resolu- 
tion as was emphasized by Veillette et ali^ The theory 
suggests though that experiments with higher momen- 



tum resolution are likely to reveal considerably sharper 
low energy spin- wave peaks. 

With regard to the shape of the broad scattering at 
higher energies, we are inclined toward a picture of scat- 
tering by quasi-one-dimensional spinons. With a bare 
interchain coupling of 0.128 meV, which we know to be 
considerably renormalized downward by quantum fluctu- 
ations, one might expect that, at energies above about 
0.3 meV, the system would behave like a collection of un- 
coupled chains. In this sense we would expect the high- 
energy response of the system to be better described by 
one-dimensional spinons than by a low-order spin-wave 
perturbation calculation such as the one that we have 
studied. 
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APPENDIX A: 1/(2S) 2 CORRECTION TO THE 
GROUND STATE ENERGY AND ORDERING 
WAVE VECTOR 

The ground state energy correction at the order 
1/(2S 1 ) 2 , denoted here as Eq\q), consists of two parts, 



e£\q)=£ 3 (Q)+£ 4 (Q). 



(Al) 



£ 4 {Q) is obtained simply by substituting the Bogoli- 
ubov transformation Eq. 1)22(1 into Ti^ and perform- 
ing the normal ordering in the terms containing two 
pairs of quasiparticle operators c and c^. In this Ap- 
pendix, we present the results of calculations in which 
the Dzyaloshinskii-Moria interaction D = 0, so that 




w 2J^-p - B >-p) — ^ — 

k,p 



(A2) 



£ 3 {Q) ensues due to the key role of the cubic inter- 
actions <|19fl in the total Hamiltonian. It is technically 
advantageous first to express (|19f) in terms of the quasi- 
particle c-operators, so that 



1,2,3 



.)(. 



J J J 



+T 2 (1,2; 3) (^^-cicW)}, (A3) 
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ri(l,2,3) = (Ci + C 2 ) (u lU2 v 3 + v ± v 2 u 3 ) 

+ (Cx + C 3 ) (M1U3W2 + «lW3M2) 

+ (C 2 + C 3 ) (u 2 U3«i + , (A4) 

r 2 (l,2;3) = (Ci + C 2 )(uiu 2 M3 + U1W2W3) 

+ (Ci - C 3 )(wiM 2 U3 + UlV 2 V 3 ) 



{C 2 — C 3 )(u±v 2 u 3 +vxu 2 v 3 ). (A5) 



At T = 0, the only contribution to £d,{Q) is given by the 
'sunrise' diagram containing two vertices Ti(q, p, — q — 
p). Using the formulas l)23|) . we arrive at the following 
expression for £3: 



1 .y^PiC^'P'-q-P)]' 



167V 



—y 



(£p + £q + £p+q) 



c 2 



£p+q 



+ C p Cq 



(A p+q + S p+q )(Ap + Bp)(A q + B q ) Ap +q -Bp +q Ap + Bp Aj + Sq 



-p+q£q£p 



Numerical calculations show that for the interaction 
constants given by Eqs. J3J), 



£ 3 {Qo)/J = -0.0425, 
£4Qo)/J= 0.0103, 



(A7) 
(A8) 



which are quite small. The 1/(2S) 2 correction to Qo can 
be computed by substituting the expansion 



Q = Qo 



2S (25)2 
into the series for the ground state energy 

Eg(Q) = s 2 j q + se£\q) + E%\Q), 



(A9) 



(A10) 



and expanding its minimum in powers of 1/(25). As a 
result, 



AQ 



(2) 



\d 2 .j Q 


-1 r 


[ dQ 2 _ 





l d 3 J Q 
2 dQ 3 



• (AQ«) 2 



dQ 2 



+2 



(AQ( 1 )) + 4< ) ^ 



dQ 



(All) 



The calculations of derivatives over Q in expressions en- 
tering Eqs. l|A6|l . (|A^|) are lengthy but straight- 
forward. Omitting the corresponding cumbersome final 
expressions, we state here that as a result of numerical 
integration: 



AQ (2) /2vr « -0.011. 



(A12) 



The value of the ground state energy at its minimum 
(S = 1/2) 



S 2 J Q + SE£ ) (Q) + E i2) (Q) 



l d 2 J Q 
8 dQ 2 



(AQW) 2 



-0.459. 



(A13) 



Note that the last term in Eq. (|A13|) is always negative and arises solely due to the shift in the ordering wave vector 

AQW. Numerically however, this term is tiny, even compared to the second order correction Eq\Qo) that itself 
gives a small contrbution to the total energy Eq. (|A13() . 
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